#include "stdafx.h"
#include "CubicSplain.h"

#pragma region "Constructors\Destructors"
CubicSplain::CubicSplain() : splines(NULL)
{

}

CubicSplain::~CubicSplain()
{
	free_mem();
}

void CubicSplain::free_mem()
{
	delete[] splines;
	splines = NULL;
}
#pragma endregion

void CubicSplain::build_spline( const double *x, const double *y, double leftGU, double rightGU, std::size_t n )
{
	free_mem();
	this->numberGridPoints = n;

	splines = new SplineTuple[n+1];
	for (std::size_t i = 0; i < n+1; ++i)
	{
		splines[i].x = x[i];
		splines[i].a = y[i];
	}
	splines[0].c = leftGU;
	splines[n].c = rightGU;

	double *alpha = new double[n];
	double *beta = new double[n];
	double A, B, C, F, h_i, h_i1, z;
	alpha[0] = 0.;
	beta[0] = leftGU;

	for (std::size_t i = 1; i < n; ++i)
	{
		h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i];
		A = h_i;
		C = (-1) * 2. * (h_i + h_i1);
		B = h_i1;
		F = (-1) * 6. * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) / h_i);
		z = (C - A * alpha[i - 1]);
		alpha[i] = B / z;
		beta[i] = (F + A * beta[i - 1]) / z;
	}

	for (std::size_t i = n - 1; i > 0; --i)
		splines[i].c = alpha[i] * splines[i + 1].c + beta[i];

	delete[] beta;
	delete[] alpha;

	for (std::size_t i = 1; i<n+1; i++)
	{
		double h_i = x[i] - x[i - 1];
		splines[i].d = (splines[i].c - splines[i - 1].c) / h_i;
		splines[i].b = h_i * (2. * splines[i].c + splines[i - 1].c) / 6. + (y[i] - y[i - 1]) / h_i;
	}
}

SplineTuple* CubicSplain::getSplineForX(double x) const
{
	SplineTuple *s;
	if (x <= splines[1].x) 
		s = &splines[1];
	else if (x >= splines[numberGridPoints-1].x) 
		s = &splines[numberGridPoints];
	else 
	{
		std::size_t first = 2, last = numberGridPoints+1;
		size_t mid;
		while (first < last)
		{
			mid = first + (last - first) / 2;
			if (x <= splines[mid].x)
				last = mid;
			else
				first = mid + 1;
		}
		s = &splines[last];
	}
	return s;
}

double CubicSplain::getF(double x) const
{
	if (!splines)
		return std::numeric_limits<double>::quiet_NaN();

	SplineTuple* s = getSplineForX(x);
	double dx = (x - s->x);
	return s->a + (s->b + (s->c / 2. + s->d * dx / 6.) * dx) * dx; 
}

double CubicSplain::getDiffOne(double x) const
{
	if (!splines)
		return std::numeric_limits<double>::quiet_NaN();
	
	SplineTuple* s = getSplineForX(x);
	double dx = (x - s->x);
	return s->b + (s->c + s->d * dx / 2.) *dx;
}

double CubicSplain::getDiffTwo(double x) const
{
	if (!splines)
		return std::numeric_limits<double>::quiet_NaN();

	SplineTuple* s = getSplineForX(x);
	double dx = (x - s->x);
	return s->c + s->d * dx;
}

void CubicSplain::getSplineCoefficients( int index, double& a, double& b, double& c, double& d, double& x )
{
	a = splines[index].a;
	b = splines[index].b;
	c = splines[index].c;
	d = splines[index].d;
	x = splines[index].x;
}
